p 


FINAL  TECHNICAL  REPORT 


AD-A237  563 


Submitted  by 


r 


D jTT?.i3r.T! o r '  ?rXr£Mrrrr  k  ; 

Kzpzz"*-r.  :cr  p-_ri ;c 

j’ • -*-d 


Prof.  Karol  Z.  Korczak,  P I. 
Mechanical  and  Aerospace  Engineering 
Case  Western  Reserve  University 
10900  Euclid  Avenue,  Cleveland,  Ohio  44106 


ONR  Program  Manager  -  James  A.  Fein 
Contract  Scientific  Officer  -  V.J.  Monacella 


1 


ABSTRACT 


The  research  project  has  been  initiated  to  develop  a  numerical  algorithm  and  conduct  a 
broad  investigation  of  flow  structures  and  their  interactions  with  moving  and  ejecting  objects.  The 
numerical  method  is  based  on  a  spectral  element  approach  applied  to  the  unsteady,  incompressible 
Navier-Stokes  equations.  The  equations  are  adopted  for  a  moving  mesh  system  that  allows  the 
mesh  to  adjust  dynamically  to  changing  geometry.  The  method  has  been  used  to  study  flow  pat¬ 
terns  and  forces  on  moving  objects.  Among  the  investigated  cases  are:  two-  and  three-dimensional 
piston  pumps,  flow  through  an  oscillating  channel,  an  oscillating  cylinder  in  a  still  fluid,  an  oscil¬ 
lating  cylinder  in  a  uniform  flow,  and  an  ejecting  object  from  a  silo. 
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INTRODUCTION 


Real  flow  cases  with  objects  deforming  or  moving  within  the  flow  are  quite  common.  Ma¬ 
chines  with  moving  parts,  structures  deforming  or  vibrating  in  the  wind  or  water  current,  artificial 
heart  pumping  blood,  opening  and  closing  valves,  aircraft  dropping  bombs,  submarines  ejecting 
torpedoes  or  launching  missiles.  These  examples  represent  but  few.  Numerical  simulations  of  such 
cases  are  not  common  due  to  relatively  high  complexity  of  the  geometry  involved,  numerical  dif¬ 
ficulties  with  handling  nonstationary  meshes,  high  complexity  of  the  unsteady  flow  fields,  and 
large  CPU  demands.  Even  experimental  investigations  are  relatively  rare  and  the  evidence  in  liter¬ 
ature  is  even  more  limited  due  to  the  fact  that  some  of  the  research  involved  classified  data  or  re¬ 
sults.  In  many  cases,  the  scientific  rewards  for  the  unavoidable  significant  research  effort  might  not 
be  as  attractive  as  in  more  traditional  areas. 

In  this  study,  the  primary  focus  was  on  phenomena  associated  with  ejection  of  objects  into  sur¬ 
rounding  fluid.  However,  several  other  cases  were  investigated  for  comparison  of  numerical  and 
experimental  results  and  to  determine  overall  possibility  of  such  investigations. 

The  spectral  element  method  [1-9]  was  adopted  for  this  investigation.  The  method  was  refor¬ 
mulated  to  allow  the  deformations  and  motion  of  the  numerical  mesh  which  followed  the  moving 
objects.  Initially,  a  piston  pump  was  simulated  (two-  and  three-dimensional)  to  provide  a  test  case 
for  the  method  and  to  allow  qualitative  comparisons  with  existing  experimental  evidence  [10,1 1]. 

The  problems  of  underwater  motion,  ejection,  and  separation  has  been  the  focus  of  several  in¬ 
vestigations  [12-20],  primarily  associated  with  military  applications  such  as  launches  from  subma¬ 
rines  of  missiles  or  torpedoes.  More  general  applications  are  associated  with  oscillating  slim 
underwater  structures  [21-29]  such  as  cylindrical  supports  of  oil  platforms  and  bridges.  That  cate¬ 
gory  includes  also  toll  buildings  and  chimneys  subjected  to  the  forces  resulting  from  air  flow. 

This  investigation  simulated  and  analyzed  many  of  those  type  of  flows  with  results  described 
in  the  following  sections.  The  analysis  in  most  cases  could  not  be  extensive  due  to  the  very  large 
CPU  time  needed  to  advance  the  simulation.  Unfortunately,  this  is  the  primaiy  limiting  factor  in 
advancing  the  numerical  analysis  for  cases  like  the  ones  mentioned  above  beyond  research  interest 
to  result  in  a  design  tool. 
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NUMERICAL  APPROACH 


The  numerical  method  implemented  for  this  project  is  based  on  the  spectral  element  method 
which  is  a  blend  of  a  spectral  approach  [1,2,3]  and  finite  element  technique.  The  spectral  element 
technique  underwent  many  modifications  from  its  initial  formulation  [4,5]  to  expand  the  realm  of 
applications.  The  approach  resembles  the  h-version  of  finite  elements  [6];  however,  due  to  spectral 
local  expansions,  it  is  much  more  accurate  and  suitable  for  complex  flow  simulations  [7-8]. 

As  in  finite  elements,  the  approach  breaks  the  domain  into  a  series  of  subdomains  mapped  from 
the  physical  space  (x,y,z)  into  local  space  (r,s,t)  where  they  are  a  regular  cubes  ranging  from  - 1  to 
+1  along  each  axis.  The  geometry,  velocity,  and  pressure  are  represented  as  tensor-product  high 
order  Lagrangian  interpolants  (based  on  Legendre  polynomials)  through  Gauss-Lobatro  points. 
Those  points  are  the  roots  of  the  first  derivative  of  the  Legendre  polynomials  and  the  end  points. 
The  equations  are  formed  for  each  subdomain  (element)  and  then  combined  into  a  global  matrix 
system  using  the  direct  stiffness  method  (as  in  finite  elements). 

The  approach  uses  the  unsteady,  incompressible  Navier-Stokes  equations  [9]  which  together 
with  mass  conservation  have  the  following  form: 


and 


5v;  3v; 

Tt 


(l-i) 


V*  V  =  0  (1.2) 

where: 

V  =  vji  +•  vn j  vyk  =  velocity 
p  ~  static  pressure 
p  =  density 

v  =  kinematic  viscosity 
Fi  ~  body  force 


-  a  - 


To  solve  the  system,  a  splitting  scheme  is  implemented  with  the  following  steps: 
-  nonlinear  step: 


or'-o 

At 


(13) 


-  pressure  step: 


K-  -V;  ) 


At 


(1.4) 


viscous  step: 


dr 


(1-5) 


The  non-linear  step  is  solved  explicitly  using  the  third-order  Adams-Bashforth  approach  (no 
boundary  conditions  are  needed  here): 


(v; 


a  i 


(1-6) 


where  Cq  =  23/32,  Ci=-16/12,  C2=5/12. 

y  /\/ 

The  system  above  perform*  satisfactory  if  the  Courm:  cumber  i— - )  >s  nor  jreaici  man  0 

The  pressure  step  is  solved  implicitly.  The  icnown  pressure  is  imposed  at  input  as  Dirichlet  bound¬ 
ary  condition  with  natural  boundary  conditions  everywhere  else.  The  viscous  step  is  solved  using 
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Crank-Nicholson  approach  with  the  required  boundary  conditions  on  velocity.  At  the  outflow,  the 
natural  boundary  conditions  are  imposed. 

To  allow  dynamic  changes  of  the  mesh  during  simulations,  the  momentum  equations  must  be 
modified  to  account  for  that  additional  motion.  To  achieve  that,  the  new  positions  of  the  deforming 
mesh  (moving  mesh  points)  must  be  known.  Together  with  the  previous  positions,  the  information 
is  sufficient  to  find  the  new  velocity  and  pressure  fields  at  mesh  points  in  the  new  locations.  The 
reformulation  of  the  momentum  equation  is  fairly  straightforward  and  does  not  add  significantly 
to  the  complexity  of  the  regular  spectral  element  codes. 

Let  x  =  f(x,  t)  to  be  the  old  position  of  the  mesh  and  x  =  f(x,  t)  to  be  the  new  position  of 
the  mesh.  The  old  nodal  velocity  would  then  be  v  =  v  (x, :)  and  the  new  nodal  velocity  would 
then  be  v  =  v  (x,  t) .  With  these,  the  momentum  equation  at  the  new  position  with  respect  to  the 
old  position  is  given  as: 


crv-  fdx  dx:  \dvt 


(1.8) 


Let  a.  = 

j 

mesh.  Withe 


UJL  ■  '*'*’*/ 

be  the  velocity  of  the  mesh  nodes  and  jL  =  be  the  local  deformation  of  the 
these,  the  final  form  of  the  momentum  equations  is  as  follows: 


3v; 

Tt 


5V; 


v^-Vs; = 


ay 


-  -r  |3?y — +  F 


(1.9) 


During  flow  simulation,  the  new  position  of  the  mesh  is  determined,  the  above  coefficients  are 
calculated  and  the  new  flow  field  is  computed.  The  extra  work  associated  with  the  moving  mesh  is 
relatively  small,  probably  no  more  than  5%  in  CPU  time.  The  benefit  of  this  approach  is  that  no 
additional  approximations  arc  necessary  and  that  the  flow  field  “follows”  the  mesh. 

One  must  be  careful  not  to  allow  the  mesh  to  deform  too  much.  If  the  deformations  are  too  ex¬ 
cessive,  the  errors  due  tc  rh^e  deformation  might  cause  serious  problems  for  accuracy  and  most 
of  ail  for  convergence  of  iterative  solvers. 
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NUMERICAL  SIMULATIONS 


The  study  involved  several  cases  of  flows  associated  with  moving  geometries.  Most  of  the  cas¬ 
es  were  investigated  to  determine  the  method’s  performance  and  ability  to  handle  such  flows.  The 
numerical  results  were  compared  when  possible  with  existing  experimental  evidence.  Both  two- 
and  three-dimensional  cases  were  considered. 


T wo-Dimensional  Piston  Pump 

The  two-dimensional  piston  pump  (Figure  1)  has  a  mesh  with  26  spectral  elements,  25  of  them 
constitute  the  inside  of  the  cylinder  and  the  26th  element  resides  in  the  inflow  neck.  Each  element 
has  125  nodes;  the  top  and  the  bottom  of  the  geometry  are  curvilinear  to  resemble  an  actual  pump. 

During  a  pumping  cycle,  the  volume  increases  2.6  times  with  respect  to  the  initial  volume.  The 
numerical  conditions  on  pressure  impose  a  fixed  pressure  at  inflow;  therefore,  an  area  around  the 
center  of  the  piston  surface  has  been  designated  as  “numerical”  inflow  with  velocity  equal  to  the 
velocity  of  the  piston  itself.  This  arrangement  fixes  an  area  with  the  reference  pressure  (here  set  to 
zero  for  convenience)  and  preventing  any  flow  through  the  piston  surface.  The  flow  may  enter  or 
exit  only  through  the  inflow  neck. 

The  motion  of  the  piston  is  enforced  to  follow  a  cosine  curve.  The  peaks  of  the  piston’s  dis¬ 
placement  correspond  to  zero  velocity.  The  maximum  velocity  occurs  in  the  midpoint  of  the  piston 
cycle.  One  piston  cycle  takes  1500  time  steps  which  corresponds  to  46.75  time  units  for  the  code. 
In  SUN  3/50  CPU  time,  it  requires  84  hours  to  run  a  full  cycle.  The  velocities  on  the  moving  piston 
surface  are  calculated  every  step  and  imposed  appropriately  for  the  flow  calculations. 

As  expected,  the  rapid  expansions  and  contractions  of  the  pump  volume  are  creaung  complex 
vortex  structures  (Figure  2).  To  avoid  the  influence  of  the  initial  conditions,  the  results  are  recorded 
starting  with  the  third  cycle  of  the  simulations. 


Three-Dimensional  Piston  Pump 

The  mesh  for  the  pump  contains  57  spectral  elements  (Figure  3)  with  125  nodes  each.  As  in  the 
two-dimensional  case,  the  mesh  has  four  layers  of  elements  with  thirteen  elements  per  layer.  The 
surface  of  the  piston  and  the  head  of  the  cylinder  chamber  are  curvilinear.  The  inlet  has  five  ele¬ 
ments  in  it.  The  pump’s  volume  expands  twice  with  respect  to  its  initial  state.  The  piston  motion 
follows  a  cosine  curve.  One  full  cycle  requires  1700  time  steps  which  corresponds  to  49.61  time 
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units.  In  SUN  4/1 10  CPU  time  units,  it  takes  600  hours  per  cycle.  Due  to  the  excessive  simulation 
time,  the  experiment  was  curried  only  through  one  and  a  half  cycle. 

The  flow  field  is  documented  in  Figure  4.  Due  to  the  acceleration  and  deceleration  of  the  piston, 
the  flow  patterns  exhibit  complex  vortex  structures  forming  in  the  cylinder  chamber. 

The  results  of  this  simulation  have  been  compared  qualitatively  with  the  experimental  results 
obtained  by  Durst,  Ernst,  and  Pareira  [10]  using  laser  doppler  anemometer.  The  experimental  re¬ 
sults  were  processed  numerically  for  visualization  [II]  and  are  almost  identical  to  those  shown  in 
Figure  5. 

The  3-D  numerical  simulation  exhibits  all  the  details  of  the  flow  demonstrated  in  the  experi¬ 
mental  investigation.  The  2-D  simulations  show  more  extensive  vortex  structure  field  due  to  the 
constraint  coming  from  the  two-dimensionality. 


The  Oscillating  Channel 

The  dynamic  motion  of  geometries  in  a  flow  results  in  variable  forces  and  moments  acting  on 
those  objects.  Those  forces  and  moments  were  recorded  in  the  simulation  of  a  flow  in  a  three-di¬ 
mensional,  rectangular,  oscillating  channel  (Figure  6). 

The  channel  consists  of  six  elements  with  343  nodes  per  element.  The  motion  imposed  on  the 
channel  is  a  sinusoidal  oscillation  of  the  exit  section  with  an  amplitude  of  one-fifth  the  channel 
height  in  a  direction  perpendicular  to  the  flow.  Forces  are  recorded  on  the  upper  and  lower  walls. 
The  moment  is  calculated  at  a  point  which  acts  as  a  hinge  for  the  oscillation  located  at  the  center 
of  the  bottom  wall  in  the  spanwise  direction.  With  a  time  step  of  0.0424,  forces  and  moments  are 
recorded  every  three  time  steps.  A  full  cycle  consists  of  300  time  steps  which  corresponds  to  nine 
CPU  hours  of  SUN  3/50. 

As  the  channel  oscillates,  the  motion  of  the  wails  imposes  a  velocity  on  fluid  near  the  wall  (Fig¬ 
ure  7).  During  the  oscillation,  the  inertia  of  the  fluid  causes  pockets  of  high  pressure  to  appear  near 
the  inlet  and  pockets  of  low  pressure  to  occur  near  the  outflow.  As  expected,  moments  and  forces 
display  sinusoidal  oscillations  synchronized  with  the  geometrical  motion  (Figure  8). 

The  results  from  this  simulation  allow  a  qualitative  evaluation  of  the  method’s  capability  to  de¬ 
termine  forces  and  moments  resulting  from  geometry-flow  structures  interactions.  For  quantitative 
evaluations  (calibration)  more  extensive  simulations  and  actual  experimental  measurements  would 
have  to  be  conducted. 
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Analysis  of  Cylinders  Oscillating  in  a  Still  Fluid  and  in  a  Uniform  Flow 


Introduction 

Analysis  of  the  flow  around  an  oscillating  cylinder  in  a  still  fluid  and  in  a  uniform  flow  shows 
that  these  cases  are  quite  different  than  flow  in  the  classical  case  of  a  cylinder  in  a  steady,  uniform 
flow.  An  oscillating  cylinder  reverses  direction,  causing  the  vortex  wake  of  one  half  cycle  to  sweep 
past  the  cylinder  during  a  subsequent  half  cycle.  This  reversal  has  a  major  effect  on  the  magnitudes 
of  fluid  induced  forces  and  also  on  the  fundamental  frequency  of  the  drag  force. 

The  interest  in  studying  flow  around  oscilladng  cylinders  has  arisen  from  many  modem  day 
construction  problems.  Vortex  induced  vibrations  are  often  the  cause  of  costly  construction  delays 
and  failures  in  offshore  structures  such  as  underwater  piping  and  cylindrical  structure  components. 
Civil  engineers  have  encountered  similar  problems  with  cylindrical  columnar  buildings  and  chim¬ 
ney  stacks  in  a  crosswind.  Force  oscillations  are  a  significant  factor  when  considering  the  fatigue 
life  of  a  structure.  Furthermore,  there  is  a  considerable  fundamental  interest  in  these  flows  and  in 
the  ways  in  which  vortices  are  generated  and  shed. 

For  years,  designers  of  offshore  platforms  typically  have  used  what  is  known  as  the  Morrison 
force  equation  [21]  to  calculate  the  wave  on  offshore  cylindrical  structures.  The  Morrison  descrip¬ 
tion  considers  the  force  to  be  composed  of  a  drag  force  (described  by  the  standard  semi-empirical 
drag  coefficient  relation)  and  an  inertia  force  (described  in  terms  of  the  inertial  coefficient  which 
is  related  to  the  hydrodynamic  added  mass  coefficient)  [22]. 

For  a  circular  cylinder,  the  in  line  force  per  unit  length  is: 

F  =  IpDC.J/M  + 

where:  p  -  fluid  density 

D  -  body  diameter 
U  -  fluid  velocity 
Cjn  -  coefficient  of  mass 
Q  -  coefficient  of  drag 

The  coefficients  and  Q  have  been  evaluated  experimentally  from  full-scale  and  laboratory 
studies  by  a  considerable  number  of  investigators.  This  is  usually  done  with  the  method  of  least 
squares  or  using  the  Fourier  averaging  over  a  cylinder  cycle  [23]. 

Keuiegan  and  Carpenter  [24]  have  found  that  the  average  values  of  and  Qj  are  functions  of 
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the  oscillating  cylinder’s  amplitude  and  diameter.  When  the  relative  flow  past  a  cylinder  is  under¬ 
going  sinusoidal  oscillations,  the  Keulegan-Carpenter  number  KC  is  defined  as  (27ta)/d;  where  a  is 
the  amplitude  of  the  relative  motion  and  d  is  the  cylinder  diameter.  Sarpkaya  [23]  has  found  that 
Qn  and  Cd  depend  not  only  on  KC,  but  also  on  a  parameter  (3  which  is  proportional  to  the  Reynolds 

number  Re:  3  =(Re)/(KC)  =  D^v,  where  T  is  the  period  of  oscillation  and  v  is  the  kinematic  vis¬ 
cosity.  This  number  is  also  known  as  the  Stokes  number,  named  after  Stokes  who  in  1851  have 
shown  that  the  force  acting  on  a  cylinder  or  a  sphere  oscillating  sinusoidally  in  a  viscous  fluid  is 
dependent  on  both  the  amplitude  of  the  oscillations  and  the  Reynolds  number. 

Williamson  [25]  has  performed  investigations  on  the  relation  between  motions  of  vortices  and 
forces  exerted  on  a  circular  cylinder  for  various  KC  numbers.  He  did  not  examine  the  three-dimen¬ 
sional  structure  of  the  flow  along  the  cylinder  axis.  Tatsumo  and  Bearman  [26]  examined  this 
three-dimensional  structure  and  concluded  that  at  low  KC  numbers  and  at  low  Stokes  numbers  the 
flow  around  cylinders  is  two-dimensional.  For  this  reason,  a  two-dimensional  numerical  simula¬ 
tions  has  a  legitimate  justification. 

Numerical  simulations  have  been  used  quite  frequently  to  calculate  flow  around  cylinders. 
Most  of  the  time,  the  cylinders  have  been  stationary  due  to  difficulties  associated  with  the  moving 
geometry.  The  major  problem  has  been  with  describing  boundary  conditions  at  me  continuously 
accelerating  or  decelerating  solid  wall  using  the  finite  element  or  finite  difference  traditional  grid 
system  [27].  The  problem  has  been  dealt  with  in  different  ways  by  many  investigators  [27,28]. 


The  Two-Dimensional  Oscillating  Cylinder 

As  indicated  above,  a  two-dimensional  flow  exists  for  small  KC  numbers.  Furthermore,  two- 
dimensional  simulations  can  be  conducted  within  a  reasonaUe  computer  time,  making  this  a  useful 
tool  for  method  verifications  and  preliminary  studies. 

In  the  case  of  this  study,  the  two-dimensional  cylinder  mesh  contains  40  spectral  elements  with 
125  nodes  per  element  (Figure  9).  The  cylinder  has  a  unit  radius  and  the  mesh  extends  radially 
away  from  the  cylinder  to  a  distance  of  eight.  This  geometric  domain  has  been  deemed  large 
enough  by  using  experimental  visualization  results.  All  secondary  and  primary  flow  patterns  pro¬ 
duced  by  the  oscillating  cylinder  will  occur  within  this  distance.  Primary  flows  are  defined  as  those 
which  occur  outside  of  the  cylinder’s  boundary  layer.  Secondary  flows  are  those  which  occur  in¬ 
side  the  cylinder’s  boundary  layer. 

A  coordinate  system  has  its  origin  in  the  center  of  the  cylinder  at  initial  position.  The  cylinder 
oscillates  along  the  y-axis.  Boundaries  parallel  to  the  y-axis  are  defined  as  inflows  and  the  bound- 
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aries  parallel  to  the  x-axis  are  prescribed  as  outflows.  Pressure  as  weil  as  velocity  are  set  to  be  zero 
at  inflows. 

Three  test  cases  have  been  chosen,  each  with  a  different  amplitude:  0.8*radius,  1.0*radius, 
2.0*radius.  According  to  the  findings  of  Sorokodum  [29],  these  three  cases  will  result  in  three  dif¬ 
ferent  flow  fields.  Parameters  for  the  test  cases  are  summarized  in  the  table  below: 


Test  #1  Test  #2  Test  #3 

Amplitude  0.8*radius  1.0*radius  2.0*radius 

KC  Number  2.5  3.2  6.3 

Stokes#  62.7  30.8  30.8 

P.eynolds  #  157.6  96.8  193.5 


In  all  three  cases  the  period  of  oscillations  is  300  time  steps.  A  sinusoidal  motion  is  imposed 
for  the  cylinder.  To  avoid  picking  up  flow  patterns  due  to  the  start  up  of  morion,  results  are  record¬ 
ed  with  the  beginning  of  the  second  oscillation  cycle.  In  SUN  3/50  CPU  time,  these  computations 
averaged  168  hours. 

The  results  of  the  simulations  are  presented  on  Figures  10,  11,  12,  and  13.  The  results  are  in 
qualitative  agreement  with  the  works  of  Williamson,  Tatsumo  and  Bearman,  and  Sarapkaya. 


The  Three-Dimensional  Oscillating  Cylinder 

The  three-dimensional  cylinder  mesh  contains  100  spectral  elements  with  125  nodes  per  ele¬ 
ment  (Tugure  14).  The  mesh  is  an  extension  of  the  two-dimensional  case  in  the  third  dimension. 
The  cylinder  is  set  to  oscillate  with  an  amplitude  equal  to  its  radius  and  period  of  oscillations  equal 
to  3.1  time  units  (300  time  steps). 

The  three-dimensional  case,  due  to  extensive  requirements  for  CPU  time,  has  been  curried  only 
for  a  limited  time  to  determine  dynamic  forces  acting  on  the  cylinder.  Two  runs  were  performed  at 
a  free  stream  Reynolds  number  of  80  and  120,  and  a  KC  number  of  3. 14.  Forces  on  the  entire  cyl¬ 
inder  were  recorded  every  two  time  steps.  The  moment  of  forces  is  calculated  with  respect  to  the 
center  of  the  cylinder.  The  results  are  shown  on  Figure  15. 

The  results  conform  to  experimental  results.  The  forces  are  plotted  as  a  function  of  time  and  as 
a  function  of  velocity.  The  relation  between  forces  and  velocities  in  that  type  of  motion  was  used 
by  Hamann  and  Daito  [22]  to  try  to  find  a  nondimensional  parameter  to  correlate  forces  due  to 
wave  action.  Their  findings  show  that  the  measured  forces  are  more  dependent  on  the  velocity  than 
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on  the  acceleration.  They  concluded  that  for  a  constant  value  of  A Jd  the  maximum  force  per  unit 
area  is  a  function  of  the  actual  velocity  and  that  the  force  is  increasing  with  the  increase  of  the  A / 
d  ratio.  The  slope  of  the  curve  from  the  results  of  the  current  investigation  is  in  agreement  with 
those  from  experiments. 


The  Separating  Object  Case 


Introduction 

An  initial  objective  to  study  full  separation  process  of  an  object  ejected  into  a  flow  field  has 
been  modified  following  suggestions  from  the  David  Taylor  Research  Center  and  focussed  on  the 
interactions  dunng  the  initial  stages  of  an  object  plowing  into  bay  and  then  into  the  surrounding 
flow  field.  The  geometry  adopted  for  this  investigation  is  shown  below. 


There  is  a  scarce  information  in  literature  on  this  subject.  Most  of  the  work  done  has  been  ex¬ 
perimental  and  theoretical  and  was  directed  towards  underwater  bodies.  The  available  papers 
present  an  evidence  that  can  be  divided  into  four  groups  of  interests:  a)  equations  of  motion  of  sub- 
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merged  bodies  [14,15],  b)  water  entry  and  exit  problems  [16,17],  c)  forces  acting  on  submerged 
bodies  [18,19]  and  d)  journal  articles  related  to  missiles  and  submarines  [13,20],  Most  of  the  re¬ 
search  is  considered  classified  and,  therefore,  with  limited  access. 

The  case  under  investigation  in  this  research  involves  the  geometry  configuration  shown  earli¬ 
er.  The  object,  initially  in  the  silo,  is  pushed  into  the  surrounding  fluid  until  the  moment  when  sep¬ 
aration  would  occur.  The  interest  here  is  on  the  flow  structures,  forces,  and  moments  acting  on  the 
object  durng  that  initial  stage  of  ejection.  The  ejection  could  be  into  a  still  fluid  or  into  a  cross- 
flow. 

Due  to  a  relatively  complex  geometry,  an  automated  mesh  generator  routine  has  been  created 
with  five  main  parameters:  the  object  (inner)  radius,  the  silo  (outer)  radius,  the  angle  of  the  silo’s 
surface  (theta),  the  height  of  the  object,  and  the  initial  height  of  the  object  in  the  silo.  The  generated 
mesh  consists  of  193  spectral  elements  (Figure  16).  The  mesh  is  constructed  in  such  a  way  that 
when  all  elements  are  in  their  unstretched  configuration,  the  object  is  fifty  six  percent  of  its  lengths 
above  the  silo.  In  that  position,  the  object’s  length  to  diameter  ratio  is  5.5:1  which  resembles  the 
proportions  for  a  Trident  C4  missile. 

The  way  the  mesh  changes  during  the  ejection  process  affects  the  accuracy  and  speed  of  the 
calculations  due  to  the  fact  that  matrix  solvers  are  based  on  iterative  schemes  and  more  distorted 
elements  require  more  iterations  to  converge  to  a  solution.  In  addition,  it  was  discovered  that  there 
might  be  a  considerable  difficulty  to  impose  properly  the  conditions  for  non-zero  velocity  on  the 
outside  flow  if  the  surface  of  the  mesh  is  complicated. 

The  object  is  assumed  to  be  pushed  out  of  the  silo  with  the  necessary  force.  At  the  same  time, 
the  fluid  will  be  pushed  through  the  annulus  into  the  surrounding  flow  field.  The  flow  of  that  fluid 
is  prescribed  by  the  following  relation  (for  a  flow  through  a  concentric  annulus  with  a  moving  inner 
wall): 


{(a*-b‘)+VwaU) 


In(-) 

a 


x  /n  (-) 
a 


Vr  =  0 

where  a  -  radius  of  the  object 
b  -  radius  of  the  silo 
Vwaii  -  object’s  velocity 
Va  -  axial  velocity  component 
Vr  -  radial  velocity  component 
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Stationary  Object  Case 

The  first  investigated  case  was  sec  with  the  following  parameters:  object  radius=0.6,  silo  radi- 
us=1.2,  object  height=6.6,  object  height  above  silo=3.76,  theta=60  degrees.  The  objea  was  kept 
stationary  and  the  flow  in  the  annulus  was  initiated  with  maximum  velocity  at  the  inlet  to  the  an¬ 
nulus  being  ramped  from  zero  to  its  normal  value  (0.5  for  this  configuration)  over  50  time  steps. 
The  ramping  follows  a  cosine  curve  for  smooth  acceleration  and  deceleration  of  the  entering  flow. 
In  addition,  the  flow  was  allowed  to  settle  from  the  impact  of  the  ramping  for  the  next  50  time  steps. 
The  data  for  analysis  was  taken  after  100  time  steps  of  the  numerical  simulation.  Figure  17  shows 
the  flow  and  pressure  fields  which  are  as  one  might  expect. 

This  case  shewed  the  main  problem  for  simulations  of  that  type:  the  overwhelming  requirement 
for  CPU  time.  On  SUN  4/1 10  the  first  100  time  steps  of  the  simulation  took  432  CPU  hours.  With 
increasing  complexity  (more  deformed  mesh,  more  resolution),  the  required  time  grows  quite  dra¬ 
matically.  For  that  reason,  it  is  very  difficult  to  achieve  developed  states  for  the  flow  and  pressure 
fields. 

The  simulation  was  repeated  for  similar  parameters  with  theta=0  (flat  silo  surface).  The  results 
are  shown  in  Figure  18.  The  velocity  and  the  pressure  are  as  expected  in  such  cases.  This  case  was 
used  as  a  starting  point  (initial  conditions)  for  the  case  with  moving  object. 


Object  Ejecting  into  Still  Flow  Field  -  Forced  Motion 


To  allow  the  objea  to  move  in  all  directions,  the  elements  in  the  mesh  must  be  able  to  deform 
accordingly.  In  this  particular  case,  the  generated  mesh  contains  two  rings  of  elements  (Figure  19) 
that  change  their  shape  to  allow  the  object  to  move.  When  the  object  moves  in  the  z-direction,  el¬ 
ements  in  the  innermost  ring  shrink  or  expand  in  the  z-direction.  Elements  in  the  outermost  ring 
then  contract  or  expand  in  order  that  their  innermost  element  faces  match  with  the  corresponding 
elements  in  the  inner  ring  (Figures  below). 


Objea  Wall 


Stationary  Elements 
X 


It  is  assumed  in  this  case  that  the  object  can  pivot  with  respect  to  the  bottom  of  the  silo  (the 
inlet.  to  the  annulus).  If  the  object  moves  in  direction  perpendicular  to  main  motion,  the  mesh  trans¬ 
forms  accordingly  to  accommodate  the  new  position.  Rotation  of  die  object  with  respect  to  its  own 
axis  does  not  cause  any  mesh  deformations. 

Several  cases  have  been  simulated  with  the  object  moving  in  and  out  of  the  silo.  That  motion 
generates  intensive  vortex  structures  in  the  vicinity  of  the  object’s  tip.  The  addition  of  forced  mo¬ 
tion  into  x-y-directions  creates  quite  compiex  flow  fields  (Figure  20)  due  to  the  fact  that  the  object 
pushes  and  drags  a  large  volume  of  fluid  during  side  motions. 

Object  Ejecting  into  Cross-Flow  -  Forced  Motion 

Due  to  the  complexity  of  the  mesh  changes  during  an  object  ejection  cycle,  the  imposition  of  a 
cross-flow  through  proper  boundary  conditions  above  the  surface  of  the  silo  exhibited  some  diffi¬ 
culties.  In  the  current  case,  the  top  of  the  computational  domain  has  been  assigned  an  impermeable 
wall  condition  to  force  the  flow  to  go  along  the  surface  of  the  silo.  The  velocity  profile  perpendic¬ 
ular  to  the  silo  surface  was  assumed  to  be  Blasius  profile  for  a  boundary  layer  flow. 

When  the  top  of  the  computational  mesh  has  fragmented  boundary  conditions  (part  with  Di- 
richlet  type,  pan  wun  Neuman  type),  the  junctions  of  the  various  types  of  the  boundary  conditions 
did  not  functioned  properly  creating  comer  type  flows  and  other  problems.  The  mesh  should  be 
constructed  in  a  way  that  will  assure  uniform  boundary  conditions  on  the  top  section  of  the  mesh 
where  the  surrounding  fluid  interacts  with  the  flow  patterns  generated  by  the  ejeenng  body. 

The  simulation  was  initially  done  for  immobilized  object  followed  by  object  ejection  cycle. 
The  flow  and  pressure  patterns  (Figure  21)  are  complex  in  these  cases  and  the  simulations  must  be 
curried  on  for  some  time  to  diminish  the  influence  of  initial  conditions.  Unfortunately,  the  simula¬ 
tions  are  very  CPU  demanding  which  prevented  the  current  investigation  from  getting  fully  devel¬ 
oped  flows  and  conducting  more  detailed  study. 


CONCLUSIONS 


The  results  of  the  simulations  indicate  that  the  numerical  flow  analysis  for  such  cases  can  be 
accurate  and  predict  important  flow  phenomena.  For  three-dimensional  cases  combined  with  mo- 
don  of  bodies,  the  formulation  becomes  quite  complex  and  requires  significant  resolution  (large 
number  of  computational  points).  When  the  motion  is  significant,  the  mesh  can  become  quite  dis¬ 
torted  and  affect  the  accuracy  and  the  speed  with  which  a  solution  at  a  given  time  step  is  achieved, 
particularly  with  respect  to  a  pressure  step.  A  care  must  be  taken  to  prevent  excessive  distortions 
of  elements  in  meshes.  With  these  aspects  in  mind,  the  numerical  simulations  serve  as  a  very  useful 
tool  in  research  and  practice  involving  cases  as  investigated  in  this  project  The  comparison  be¬ 
tween  the  simulations  and  the  (when  available)  experimental  data  showed  good  qualitative  agree¬ 
ment. 

All  simulations  in  the  cases  above  were  for  situations  in  which  the  motions  were  forced  on  the 
moving  objects  and  the  resulting  forces  and  moments  from  the  body-fluid  interactions  calculated. 
However,  knowing  the  forces  and  the  moments  (torques)  acting  on  a  moving  body,  it  is  fairly 
straightforward  to  develop  a  dynamic  equations  of  motion  for  that  body  allows  to  determine  the 
trajectory  (path).  That  is  essential  for  unconstrained  bodies  such  as  ejected  objects.  A  preliminary 
program  was  developed  for  that  purpose  and  initially  tested.  Due  to  the  limited  data  from  ejecting 
object  simulations  and  time  constraints,  the  dynamic  motion  solver  has  not  been  completed. 

The  main  factor  limiting  the  scope  of  this  research  was  the  large  demand  for  CPU  time,  partic¬ 
ularly  for  three-dimensional  cases  involving  significant  object  motions  (mesh  distortions).  The  use 
of  supercomputers,  vector  processors,  and  multiprocessor  machines  could  improve  that  situations. 


RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

There  is  a  significant  interest  on  pan  of  the  industry  for  simulations  involving  moving  objects 
and  particularly  for  cases  that  require  the  simulation  not  only  of  the  flow  field  but  also  the  dynamic 
behavior  of  the  body  resulting  from  the  external  fluid  forces.  However,  at  present  it  seams  that  for 
realistic  and  accurate  simulations  of  such  cases,  the  required  size  of  the  computational  model  re¬ 
sults  in  very  large  computing  times.  In  addition,  cases  with  large  Reynolds  numbers  and  highly  un¬ 
steady  further  complicate  the  situation.  Therefore,  with  the  current  computational  capabilities 
(primarily  supercomputers  or  comparable  machines)  one  can  approach  only  limited  cases.  Even¬ 
tually,  when  the  barriers  of  computer  speed  will  be  overcome  by  technological  developments  and 
new  solvers  for  matrix  equations  (the  main  CPU  time  consumers)  will  be  developed,  the  successful 
research  developments  in  this  area  will  offer  scientific  rewards  as  well  as  technological  payoffs. 
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List  of  Figures 


Figure  1.  Schematic  drawing  of  the  computational  domain  for  the  two-dimensional  piston 
pump. 

Figure  2  (a-e)  Plots  of  velocity  vectors  and  pressure  for  the  two-dimensional  pump  cycle.  The 
period  of  the  cycle  is  identified  for  every  plot  at  the  bottom  left  comer.  Vector  scales  arc  located 
at  the  lower  right  comer  (scale  is  with  respect  to  units  identified  on  the  axis).  For  pressure,  open 
ended  arrows  indicate  negative  pressure,  close  ended  arrows  indicate  positive  pressure. 

Figure  3.  The  computational  domain  for  the  three-dimensional  piston  pump  displaying  all  57 
elements.  Schematic  drawing  of  the  piston  pump  is  identical  as  on  Figure  1. 

Figure  4.  (a-f)  Plots  of  velocity  vectors  and  pressure  contours  for  the  three-dimensional  piston 
pump  during  pump  cycle.  Values  between  contour  lines  are  identified  at  lower  right  comer  of  pres¬ 
sure  plots.  (The  first  number  corresponds  to  the  difference  between  negative  pressure  lines  (dotted) 
and  the  second  number  corresponds  to  the  difference  between  positive  pressure  lines  (solid  lines). 

Figure  5.  Vorncity  plots  from  an  experiment  using  laser  doppler  anemometer.  Sine  curves  at 
the  top  of  each  plot  identify  the  position  of  the  piston  in  the  pump’s  cycle. 

Figure  6.  Schematic  drawings  of  the  three-dimensional  oscillating  channel  geometry. 

Figure  7  (a,b).  Plots  of  velocity  vectors  and  pressure  contours  during  one  cycle  in  an  oscillat¬ 
ing  channel  flow.  Sine  curves  next  to  each  plot  represent  the  position  of  the  channel  in  the  cycle. 

Figure  8  (a,b).  Plots  of  forces  and  moments  on  the  walls  of  the  oscillating  channel.  The  x  rep¬ 
resents  the  number  of  time  steps  divided  by  three. 

Figure  9.  Schematic  drawing  of  the  two-dimensional  oscillating  cylinder. 

Figure  10(a-c).  Velocity  vectors  for  the  two-dimensional  cylinder  oscillating  at  an  amplitude 
equal  to  0.8*Tadius.  Sine  curves  at  the  top  of  each  figure  identify  the  location  of  the  cylinder  in  the 
cycle. 

Figure  ll(a-c).  Velocity  vectors  for  the  two-dimensional  cylinder  oscillating  at  an  amplitude 
equal  to  1.0*radius.  Sine  curves  at  the  top  of  each  figure  identify  the  location  of  the  cylinder  in  the 
cycle. 
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Figure  12(a-c).  Velocity  vectors  for  the  two-dimensional  cylinder  oscillating  at  an  amplitude 
equal  to  2.0*radius.  Sine  curves  at  the  top  of  each  figure  identify  the  location  of  the  cylinder  in  the 
cycle. 

Figure  13(a-c).  Pressure  for  the  two-dimensional  cylinder  oscillating  at  an  amplitude  equal 
to  i.0*radius  (opened  arrows  indicate  negative  pressure,  closed  arrows  indicate  positive  pressure). 
Sine  curves  at  the  top  of  each  figure  identify  the  location  of  the  cylinder  in  the  cycle. 

Figure  14.  Computational  domain  for  the  three-dimensional  oscillating  cylinder.  Only  ele¬ 
ment  edges  are  plotted  for  clarity. 

Figure  15.  Forces  and  moments  recorded  on  the  three-dimensional  oscillating  cylinder  at  a 
Reynolds  number  of  120. 

Figure  16.  Schematic  drawing  of  the  separating  objects  configuration. 

Figure  17.  Plots  of  velocity  vectors  and  pressure  for  the  separating  objects  case  for  theta  equal 
to  60  degrees.  (Plots  are  for  flow  after  100  time  steps  of  simulation.) 

Figure  18.  Plots  of  velocity  vectors  and  pressure  for  the  separating  objects  case  for  theta  equal 
to  zero  degrees.  For  pressure,  open  arrows  indicate  negative  pressure,  closed  arrows  indicate  pos¬ 
itive  pressure,  the  length  of  the  arrows  is  proportional  to  pressure  magnitude.  (Plots  are  for  flow 
after  100  time  steps  of  simulation.) 

Figure  19.  Example  of  the  motion  of  the  object  in  the  z-direction.  Elements  in  the  figure  on 
the  left  are  in  the  unstretched  position.  The  right  figure  shows  the  object  in  the  silo  with  elements 
stretched  to  accommodate  that  position. 

Figure  20.  Flow  field  created  by  an  object  moving  in  the  y-direction  as  shown  at  the  top  of  the 
figure.  Velocity  vectors  and  pressure  are  plotted  next  to  it  For  pressure,  open  arrows  indicate 
negative  pressure,  closed  arrows  indicate  positive  pressure,  the  length  of  the  arrows  is  proportional 
to  pressure  magnitude. 

Figure  21.  Plots  of  velocity  vectors  and  pressure  for  the  separating  objects  case  with  theta 
equal  to  10  degree.  For  pressure,  open  arrows  indicate  negative  pressure,  closed  arrows  indicate 
positive  pressure,  the  length  of  the  arrows  is  proportional  to  pressure  magnitude.  (Plots  represent 
flow  fieid  after  100  time  steps  of  simulation.) 
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Piston  Valve  Neck  (How  Inlet/Outlet) 


Figure  1 


Schematic  Drawing  of  the  Two-Dimension  Pump 
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Figure  2  (b)  Two-Dimensional  P  ump  Velocity  and  Pressure  Cycle  Plots 
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Fgure  2  (cj  Two- Dimensional  Pomp  Velocity  and  Pressure  Cycle  Plots 
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Figure  2  (d)  Two- Dimensional  Pump  Velocity  and  Pressire  Cycle  Plots 
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Figure  2  (e)  Two-Dimensional  Pump  Velocity  and  Pressure  Cycle  Plots 
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Figure  3 
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The  Three- Dimensional  Pump 
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Figure  4  (a) 


Three-Dimensional  Pump  Velocity  and  Pressure  Cycle  Plots 
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Figure  4  (c)  Three-Dimensional  Pump  Velociry  and  Pressure  Cycle  Plots 
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Figure  4  (d)  Three-Dimensional  Pump  Velocity  and  Pressure  Cycle  Plots 
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Figure  4  (c)  Three-Dimensional  Pump  Velocity  and  Pressure  Cycle  Plots 
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Figure  4  (f)  Three-Dimensional  Pump  Velocity  and  Pressure  Cycle  Plots 
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Flow  Inlet 


Flow  Outlet 
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Figure  6  The  Three- Dimensional  Oscillating  Channel 


Figure  7  (b)  Oscillating  Channel  Velocity  and  Pressure  Plots 
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Data  from  'Oscillating  Channel' 
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(Origin  at  Center  of  Cylinder) 


Figure  9  Schematic  Drawing  of  the  Two-Dimensional  Oscillating  Cylinder 
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Figtire  13  (a)  Pressure 


vectors  for  die  Oscillating  Cylinder  (Axnpi.  =  Rad.) 
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Figure  13(b)  Pressure  vectors  for  the  Oscillating  Cylinder  (AmpL  =  Rad.) 
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Figure  13(c)  Pressure  veaors  for  die  Oscillanng  Cylinder  (Ai 


pi.  =  Rad.) 
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Figure  14  The  Three- Dimensional  Oscillating  Cylinder  Mesh 
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Flow  in  an  Annulus  with  a  Moving  Inner  WalL 


Figure  16 


Schematic  Drawing  of  the  Separanng  Object  Mesh 
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Figure  13  Velocity  and  Pressure  Plots  of Theta  =  0  Degrees 
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Figure  19  Examp ie  of  Mesh  Element  Modon 
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Figure  Zi  Velocity  and  Pressure  Plots  of  Theta  =  10  Degrees 
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